---
title: "Power Analysis"
author: "Kaylyn Jackson Schiff, Daniel Schiff, and Natalia Bueno"
date: '2022'
output:
  html_document:
    df_print: paged
editor_options:
  chunk_output_type: console
---

#####Setup Chunk
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(rmarkdown)
library(tidyverse)
library(estimatr)
library(car)
options(scipen=999)
```


#####Estimates from Previous Studies
```{r previous estimates}
ld <- readRDS("data/df_clean.rds")
ld2 <- readRDS("data/df_followup_clean.rds")

mean_1_control <- mean(ld$support[ld$alleg_treatment=="Control"])
sd_1_control <- sd(ld$support[ld$alleg_treatment=="Control"])
ld_text_reg <- lm_robust(data = ld %>% filter(media_format=="Text"), 
                   support ~ alleg_treatment +
                     party + gender + race + age + education + income + region + media_literacy + digital_literacy)
te_1_iu_text <- coef(ld_text_reg)[2] %>% unname()
te_1_or_text <- coef(ld_text_reg)[3] %>% unname()
ld_video_reg <- lm_robust(data = ld %>% filter(media_format=="Video"), 
                   support ~ alleg_treatment +
                     party + gender + race + age + education + income + region + media_literacy + digital_literacy)
te_1_iu_video <- coef(ld_video_reg)[2] %>% unname()
te_1_or_video <- coef(ld_video_reg)[3] %>% unname()

mean_1_control_nodo <- mean(ld$support_nodonation[ld$alleg_treatment=="Control"])
sd_1_control_nodo <- sd(ld$support_nodonation[ld$alleg_treatment=="Control"])
ld_text_reg_nodo <- lm_robust(data = ld %>% filter(media_format=="Text"), 
                   support_nodonation ~ alleg_treatment +
                     party + gender + race + age + education + income + region + media_literacy + digital_literacy)
te_1_iu_text_nodo <- coef(ld_text_reg_nodo)[2] %>% unname()
te_1_or_text_nodo <- coef(ld_text_reg_nodo)[3] %>% unname()
ld_video_reg_nodo <- lm_robust(data = ld %>% filter(media_format=="Video"), 
                   support_nodonation ~ alleg_treatment +
                     party + gender + race + age + education + income + region + media_literacy + digital_literacy)
te_1_iu_video_nodo <- coef(ld_video_reg_nodo)[2] %>% unname()
te_1_or_video_nodo <- coef(ld_video_reg_nodo)[3] %>% unname()

mean_2_control <- 0
sd_2_control <- 1

ld_2_iu_reg <- lm_robust(data = ld2, 
                   support_exp1 ~ alleg_treatment_1 +
                     party + gender + race + age + education + income + region + media_literacy + digital_literacy)
ld_2_iu_reg_nodo <- lm_robust(data = ld2, 
                   support_exp1_nodonation ~ alleg_treatment_1 +
                     party + gender + race + age + education + income + region + media_literacy + digital_literacy)
ld_2_or_reg <- lm_robust(data = ld2, 
                   support_exp2 ~ alleg_treatment_2 +
                     party + gender + race + age + education + income + region + media_literacy + digital_literacy)
ld_2_or_reg_nodo <- lm_robust(data = ld2, 
                   support_exp2_nodonation ~ alleg_treatment_2 +
                     party + gender + race + age + education + income + region + media_literacy + digital_literacy)

library(meta)
ates <- c(summary(ld_text_reg)$coefficients[2,1], summary(ld_2_iu_reg)$coefficients[2,1]) %>% unname
ses <- c(summary(ld_text_reg)$coefficients[2,2], summary(ld_2_iu_reg)$coefficients[2,2]) %>% unname
ns <- c(nrow(ld), nrow(ld2))
support_IU <- cbind(ates, ses, ns) %>% as_tibble()
te_2_iu <- metagen(data = support_IU, TE = ates, seTE = ses, n.e = ns)$TE.fixed

ates <- c(summary(ld_text_reg)$coefficients[3,1], summary(ld_2_or_reg)$coefficients[2,1]) %>% unname
ses <- c(summary(ld_text_reg)$coefficients[3,2], summary(ld_2_or_reg)$coefficients[2,2]) %>% unname
ns <- c(nrow(ld), nrow(ld2))
support_OR <- cbind(ates, ses, ns) %>% as_tibble()
te_2_or <- metagen(data = support_OR, TE = ates, seTE = ses, n.e = ns)$TE.fixed

ates <- c(summary(ld_text_reg_nodo)$coefficients[2,1], summary(ld_2_iu_reg_nodo)$coefficients[2,1]) %>% unname
ses <- c(summary(ld_text_reg_nodo)$coefficients[2,2], summary(ld_2_iu_reg_nodo)$coefficients[2,2]) %>% unname
ns <- c(nrow(ld), nrow(ld2))
support_IU <- cbind(ates, ses, ns) %>% as_tibble()
te_2_iu_nodo <- metagen(data = support_IU, TE = ates, seTE = ses, n.e = ns)$TE.fixed

ates <- c(summary(ld_text_reg_nodo)$coefficients[3,1], summary(ld_2_or_reg_nodo)$coefficients[2,1]) %>% unname
ses <- c(summary(ld_text_reg_nodo)$coefficients[3,2], summary(ld_2_or_reg_nodo)$coefficients[2,2]) %>% unname
ns <- c(nrow(ld), nrow(ld2))
support_OR <- cbind(ates, ses, ns) %>% as_tibble()
te_2_or_nodo <- metagen(data = support_OR, TE = ates, seTE = ses, n.e = ns)$TE.fixed

ld_reg <- lm_robust(data = ld, 
                   support ~ alleg_treatment +
                     party + gender + race + age + education + income + region + media_literacy + digital_literacy)
te_1_iu <- coef(ld_reg)[2] %>% unname()
te_1_or <- coef(ld_reg)[3] %>% unname()

ld_2_reg_iu <- lm_robust(data = ld2, 
                   support_exp1 ~ alleg_treatment_1 +
                     party + gender + race + age + education + income + region + media_literacy + digital_literacy)
te_2_iu_study_2 <- coef(ld_2_reg_iu)[2] %>% unname()
ld_2_reg_or <- lm_robust(data = ld2, 
                   support_exp2 ~ alleg_treatment_2 +
                     party + gender + race + age + education + income + region + media_literacy + digital_literacy)
te_2_or_study_2 <- coef(ld_2_reg_or)[2] %>% unname()

```


#####Set Parameters for Power Analysis
```{r establish parameters for power analysis 1}
alpha <- 0.05 # Standard significance level 
possible.n <- seq(from=500, to=4000, by=100) # The sample sizes we'll be considering 
sims <- 1000 # Number of simulated treatment allocations to conduct for each N 

powers.iu <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
powers.or <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
```


####Simulate Results for Politican Response on Support Outcome (Liar's Dividend) -- Study 1 Replication, Full Support Index
```{r loop of lms to simulate results 1}
#Loop to vary the number of subjects
for (j in 1:length(possible.n)){ N <- possible.n[j] # Pick the sample size

Y0 <- rnorm(n=N, mean=mean_1_control, sd=sd_1_control) # control potential outcome 
tau1 <- te_1_iu_text # Hypothesize treatment effect for IU, text
tau2 <- te_1_or_text # Hypothesize treatment effect for OR, text
tau3 <- te_1_iu_video # Hypothesize treatment effect for IU, video
tau4 <- te_1_or_video # Hypothesize treatment effect for OR, video
Y1 <- Y0 + tau1 # treatment potential outcome for IU, text
Y2 <- Y0 + tau2 # treatment potential outcome for OR, text
Y3 <- Y0 + tau3 # treatment potential outcome for IU, video
Y4 <- Y0 + tau4 # treatment potential outcome for OR, video

significant.lm.iu <- rep(NA, sims) # Empty object to count significant experiments 
significant.lm.or <- rep(NA, sims) # Empty object to count significant experiments 

#Inner loop to conduct experiments "sims" times over for each N
for (i in 1:sims){
  
  Z.sim <-  sample(0:4, N, replace = T, prob = c(2/6, 1/6, 1/6, 1/6, 1/6)) # Do a random assignment 
  Y.sim <- case_when(Z.sim == 0 ~ Y0, 
                     Z.sim == 1 ~ Y1,
                     Z.sim == 2 ~ Y2,
                     Z.sim == 3 ~ Y3,
                     Z.sim == 4 ~ Y4) # Reveal outcomes according to assignment 

  fit.lm <- lm_robust(Y.sim ~ factor(Z.sim)) # Do analysis (Simple regression)
  
  p.lm <- summary(fit.lm)$coefficients[2,4] # Extract p-values  
  p.lm2 <- summary(fit.lm)$coefficients[3,4] # Extract p-values 
  significant.lm.iu[i] <- (p.lm <= alpha) # Determine significance according to p <= 0.05
  significant.lm.or[i] <- (p.lm2 <= alpha) # Determine significance according to p <= 0.05
  
}
  powers.iu[j] <- mean(significant.lm.iu) # store average success rate (power) for each N 
  powers.or[j] <- mean(significant.lm.or) # store average success rate (power) for each N 
}

```

####Plot Results -- Study 1 Replication, Full Support Index
```{r plot results 1}
png("Figures/power_followups_rep_iu.png")
plot(possible.n, powers.iu, type = "b", main = "Power to Detect Liar's Dividend (Text, IU) by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers.iu, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
#Need 2100

png("Figures/power_followups_rep_or.png")
plot(possible.n, powers.or, type = "b", main = "Power to Detect Liar's Dividend (Text, OR) by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers.or, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
#Need 3200
```



#####Set Parameters for Power Analysis
```{r establish parameters for power analysis 2}
alpha <- 0.05 # Standard significance level 
possible.n <- seq(from=500, to=4000, by=100) # The sample sizes we'll be considering 
sims <- 1000 # Number of simulated treatment allocations to conduct for each N 

powers.iu <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
powers.or <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
```


####Simulate Results for Politican Response on Support Outcome (Liar's Dividend) -- Study 1 Replication, Support Index without Donations
```{r loop of lms to simulate results 2}
#Loop to vary the number of subjects
for (j in 1:length(possible.n)){ N <- possible.n[j] # Pick the sample size

Y0 <- rnorm(n=N, mean=mean_1_control_nodo, sd=sd_1_control_nodo) # control potential outcome 
tau1 <- te_1_iu_text_nodo # Hypothesize treatment effect for IU, text
tau2 <- te_1_or_text_nodo # Hypothesize treatment effect for OR, text
tau3 <- te_1_iu_video_nodo # Hypothesize treatment effect for IU, video
tau4 <- te_1_or_video_nodo # Hypothesize treatment effect for OR, video
Y1 <- Y0 + tau1 # treatment potential outcome for IU, text
Y2 <- Y0 + tau2 # treatment potential outcome for OR, text
Y3 <- Y0 + tau3 # treatment potential outcome for IU, video
Y4 <- Y0 + tau4 # treatment potential outcome for OR, video

significant.lm.iu <- rep(NA, sims) # Empty object to count significant experiments 
significant.lm.or <- rep(NA, sims) # Empty object to count significant experiments 

#Inner loop to conduct experiments "sims" times over for each N
for (i in 1:sims){
  
  Z.sim <-  sample(0:2, N, replace = T, prob = c(1/3, 1/3, 1/3)) # Do a random assignment 
  Y.sim <- case_when(Z.sim == 0 ~ Y0, 
                     Z.sim == 1 ~ Y1,
                     Z.sim == 2 ~ Y2) # Reveal outcomes according to assignment 

  fit.lm <- lm_robust(Y.sim ~ factor(Z.sim)) # Do analysis (Simple regression)
  
  p.lm <- summary(fit.lm)$coefficients[2,4] # Extract p-values  
  p.lm2 <- summary(fit.lm)$coefficients[3,4] # Extract p-values 
  significant.lm.iu[i] <- (p.lm <= alpha) # Determine significance according to p <= 0.05
  significant.lm.or[i] <- (p.lm2 <= alpha) # Determine significance according to p <= 0.05
  
}
  powers.iu[j] <- mean(significant.lm.iu) # store average success rate (power) for each N 
  powers.or[j] <- mean(significant.lm.or) # store average success rate (power) for each N 
}

```

####Plot Results -- Study 1 Replication, Support Index without Donations
```{r plot results 2}
png("Figures/power_followups_rep_iu_nodo.png")
plot(possible.n, powers.iu, type = "b", main = "Power to Detect Liar's Dividend (Text, IU) by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers.iu, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
#Need 1100

png("Figures/power_followups_rep_or_nodo.png")
plot(possible.n, powers.or, type = "b", main = "Power to Detect Liar's Dividend (Text, OR) by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers.or, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
#Need 2200
```



#####Set Parameters for Power Analysis
```{r establish parameters for power analysis 3}
alpha <- 0.05 # Standard significance level 
possible.n <- seq(from=500, to=4000, by=100) # The sample sizes we'll be considering 
sims <- 1000 # Number of simulated treatment allocations to conduct for each N 

powers.iu <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
powers.or <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
```


####Simulate Results for Politican Response on Support Outcome (Liar's Dividend) -- Lucid Follow-up, Full Support Index
```{r loop of lms to simulate results 3}
#Loop to vary the number of subjects
for (j in 1:length(possible.n)){ N <- possible.n[j] # Pick the sample size

Y0 <- rnorm(n=N, mean=mean_2_control, sd=sd_2_control) # control potential outcome 
tau1 <- te_2_iu # Hypothesize treatment effect for IU, text
tau2 <- te_2_or # Hypothesize treatment effect for OR, text
Y1 <- Y0 + tau1 # treatment potential outcome for IU, text
Y2 <- Y0 + tau2 # treatment potential outcome for OR, text

significant.lm.iu <- rep(NA, sims) # Empty object to count significant experiments 
significant.lm.or <- rep(NA, sims) # Empty object to count significant experiments 

#Inner loop to conduct experiments "sims" times over for each N
for (i in 1:sims){
  
  Z.sim <-  sample(0:2, N, replace = T, prob = c(1/3, 1/3, 1/3)) # Do a random assignment 
  Y.sim <- case_when(Z.sim == 0 ~ Y0, 
                     Z.sim == 1 ~ Y1,
                     Z.sim == 2 ~ Y2) # Reveal outcomes according to assignment 

  fit.lm <- lm_robust(Y.sim ~ factor(Z.sim)) # Do analysis (Simple regression)
  
  p.lm <- summary(fit.lm)$coefficients[2,4] # Extract p-values  
  p.lm2 <- summary(fit.lm)$coefficients[3,4] # Extract p-values 
  significant.lm.iu[i] <- (p.lm <= alpha) # Determine significance according to p <= 0.05
  significant.lm.or[i] <- (p.lm2 <= alpha) # Determine significance according to p <= 0.05
  
}
  powers.iu[j] <- mean(significant.lm.iu) # store average success rate (power) for each N 
  powers.or[j] <- mean(significant.lm.or) # store average success rate (power) for each N 
}

```

####Plot Results -- Lucid Follow-up, Full Support Index
```{r plot results 3}
png("Figures/power_followups_lucid_iu.png")
plot(possible.n, powers.iu, type = "b", main = "Power to Detect Liar's Dividend (Text, IU) by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers.iu, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
#Need 4000

png("Figures/power_followups_lucid_or.png")
plot(possible.n, powers.or, type = "b", main = "Power to Detect Liar's Dividend (Text, OR) by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers.or, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
#Need 2000
```



#####Set Parameters for Power Analysis
```{r establish parameters for power analysis 4}
alpha <- 0.05 # Standard significance level 
possible.n <- seq(from=500, to=4000, by=100) # The sample sizes we'll be considering 
sims <- 1000 # Number of simulated treatment allocations to conduct for each N 

powers.iu <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
powers.or <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
```


####Simulate Results for Politican Response on Support Outcome (Liar's Dividend) -- Lucid Follow-up, Support Index without Donations
```{r loop of lms to simulate results 4}
#Loop to vary the number of subjects
for (j in 1:length(possible.n)){ N <- possible.n[j] # Pick the sample size

Y0 <- rnorm(n=N, mean=mean_2_control, sd=sd_2_control) # control potential outcome 
tau1 <- te_2_iu_nodo # Hypothesize treatment effect for IU, text
tau2 <- te_2_or_nodo # Hypothesize treatment effect for OR, text
Y1 <- Y0 + tau1 # treatment potential outcome for IU, text
Y2 <- Y0 + tau2 # treatment potential outcome for OR, text

significant.lm.iu <- rep(NA, sims) # Empty object to count significant experiments 
significant.lm.or <- rep(NA, sims) # Empty object to count significant experiments 

#Inner loop to conduct experiments "sims" times over for each N
for (i in 1:sims){
  
  Z.sim <-  sample(0:2, N, replace = T, prob = c(1/3, 1/3, 1/3)) # Do a random assignment 
  Y.sim <- case_when(Z.sim == 0 ~ Y0, 
                     Z.sim == 1 ~ Y1,
                     Z.sim == 2 ~ Y2) # Reveal outcomes according to assignment 

  fit.lm <- lm_robust(Y.sim ~ factor(Z.sim)) # Do analysis (Simple regression)
  
  p.lm <- summary(fit.lm)$coefficients[2,4] # Extract p-values  
  p.lm2 <- summary(fit.lm)$coefficients[3,4] # Extract p-values 
  significant.lm.iu[i] <- (p.lm <= alpha) # Determine significance according to p <= 0.05
  significant.lm.or[i] <- (p.lm2 <= alpha) # Determine significance according to p <= 0.05
  
}
  powers.iu[j] <- mean(significant.lm.iu) # store average success rate (power) for each N 
  powers.or[j] <- mean(significant.lm.or) # store average success rate (power) for each N 
}

```

####Plot Results -- Lucid Follow-up, Support Index without Donations
```{r plot results 4}
png("Figures/power_followups_lucid_iu_nodo.png")
plot(possible.n, powers.iu, type = "b", main = "Power to Detect Liar's Dividend (Text, IU) by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers.iu, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
#Need 2900

png("Figures/power_followups_lucid_or_nodo.png")
plot(possible.n, powers.or, type = "b", main = "Power to Detect Liar's Dividend (Text, OR) by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers.or, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
#Need 1800
```



#####Set Parameters for Power Analysis
```{r establish parameters for power analysis 5}
alpha <- 0.05 # Standard significance level 
possible.n <- seq(from=500, to=4000, by=100) # The sample sizes we'll be considering 
sims <- 1000 # Number of simulated treatment allocations to conduct for each N 

powers.iu <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
powers.or <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
```


####Simulate Results for Politican Response on Support Outcome (Liar's Dividend) -- Study 1 Replication, Full Support Index, Video vs. Text
```{r loop of lms to simulate results 5}
#Loop to vary the number of subjects
for (j in 1:length(possible.n)){ N <- possible.n[j] # Pick the sample size

Y0 <- rnorm(n=N, mean=mean_1_control, sd=sd_1_control) # control potential outcome 
tau1 <- te_1_iu_text # Hypothesize treatment effect for IU, text
tau2 <- te_1_or_text # Hypothesize treatment effect for OR, text
tau3 <- te_1_iu_video # Hypothesize treatment effect for IU, video
tau4 <- te_1_or_video # Hypothesize treatment effect for OR, video
Y1 <- Y0 + tau1 # treatment potential outcome for IU, text
Y2 <- Y0 + tau2 # treatment potential outcome for OR, text
Y3 <- Y0 + tau3 # treatment potential outcome for IU, video
Y4 <- Y0 + tau4 # treatment potential outcome for OR, video

significant.lm.iu <- rep(NA, sims) # Empty object to count significant experiments 
significant.lm.or <- rep(NA, sims) # Empty object to count significant experiments 

#Inner loop to conduct experiments "sims" times over for each N
for (i in 1:sims){
  
  Z.sim <-  sample(0:4, N, replace = T, prob = c(2/6, 1/6, 1/6, 1/6, 1/6)) # Do a random assignment 
  Y.sim <- case_when(Z.sim == 0 ~ Y0, 
                     Z.sim == 1 ~ Y1,
                     Z.sim == 2 ~ Y2,
                     Z.sim == 3 ~ Y3,
                     Z.sim == 4 ~ Y4) # Reveal outcomes according to assignment 

  fit.lm <- lm_robust(Y.sim ~ factor(Z.sim)) # Do analysis (Simple regression)
  
  p.lm <- linearHypothesis(fit.lm, "factor(Z.sim)1 = factor(Z.sim)3")[2,4] # Extract p-values  
  p.lm2 <- linearHypothesis(fit.lm, "factor(Z.sim)2 = factor(Z.sim)4")[2,4] # Extract p-values 
  significant.lm.iu[i] <- (p.lm <= alpha) # Determine significance according to p <= 0.05
  significant.lm.or[i] <- (p.lm2 <= alpha) # Determine significance according to p <= 0.05
  
}
  powers.iu[j] <- mean(significant.lm.iu) # store average success rate (power) for each N 
  powers.or[j] <- mean(significant.lm.or) # store average success rate (power) for each N 
}

```

####Plot Results -- Study 1 Replication, Full Support Index, Video vs. Text
```{r plot results 5}
png("Figures/power_followups_rep_iu_vid_text.png")
plot(possible.n, powers.iu, type = "b", main = "Deepfakes Hypothesis Power (IU) by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers.iu, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
#Need 3200

png("Figures/power_followups_rep_or_vid_text.png")
plot(possible.n, powers.or, type = "b", main = "Deepfakes Hypothesis Power (OR) by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers.or, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
#Need 2900
```



#####Set Parameters for Power Analysis
```{r establish parameters for power analysis 6}
alpha <- 0.05 # Standard significance level 
possible.n <- seq(from=3000, to=7000, by=100) # The sample sizes we'll be considering 
sims <- 1000 # Number of simulated treatment allocations to conduct for each N 

powers.iu <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
powers.or <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
```


####Simulate Results for Politican Response on Support Outcome (Liar's Dividend) -- Study 1 Replication, Full Support Index, Text and Video Pooled
```{r loop of lms to simulate results 6}
#Loop to vary the number of subjects
for (j in 1:length(possible.n)){ N <- possible.n[j] # Pick the sample size

Y0 <- rnorm(n=N, mean=mean_1_control, sd=sd_1_control) # control potential outcome 
tau1 <- te_1_iu # Hypothesize treatment effect for IU
tau2 <- te_1_or # Hypothesize treatment effect for OR
Y1 <- Y0 + tau1 # treatment potential outcome for IU
Y2 <- Y0 + tau2 # treatment potential outcome for OR

significant.lm.iu <- rep(NA, sims) # Empty object to count significant experiments 
significant.lm.or <- rep(NA, sims) # Empty object to count significant experiments 

#Inner loop to conduct experiments "sims" times over for each N
for (i in 1:sims){
  
  Z.sim <-  sample(0:2, N, replace = T, prob = c(1/3, 1/3, 1/3)) # Do a random assignment 
  Y.sim <- case_when(Z.sim == 0 ~ Y0, 
                     Z.sim == 1 ~ Y1,
                     Z.sim == 2 ~ Y2) # Reveal outcomes according to assignment 

  fit.lm <- lm_robust(Y.sim ~ factor(Z.sim)) # Do analysis (Simple regression)
  
  p.lm <- summary(fit.lm)$coefficients[2,4] # Extract p-values  
  p.lm2 <- summary(fit.lm)$coefficients[3,4] # Extract p-values 
  significant.lm.iu[i] <- (p.lm <= alpha) # Determine significance according to p <= 0.05
  significant.lm.or[i] <- (p.lm2 <= alpha) # Determine significance according to p <= 0.05
  
}
  powers.iu[j] <- mean(significant.lm.iu) # store average success rate (power) for each N 
  powers.or[j] <- mean(significant.lm.or) # store average success rate (power) for each N 
}

```

####Plot Results -- Study 1 Replication, Full Support Index, Text and Video Pooled
```{r plot results 6}
png("Figures/power_followups_rep_iu_pooled.png")
plot(possible.n, powers.iu, type = "b", main = "Power to Detect Liar's Dividend (IU) by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers.iu, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
#Need 5900

png("Figures/power_followups_rep_or_pooled.png")
plot(possible.n, powers.or, type = "b", main = "Power to Detect Liar's Dividend (OR) by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers.or, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
#Need over 7000
```



#####Set Parameters for Power Analysis
```{r establish parameters for power analysis 7}
alpha <- 0.05 # Standard significance level 
possible.n <- seq(from=2000, to=7000, by=100) # The sample sizes we'll be considering 
sims <- 500 # Number of simulated treatment allocations to conduct for each N 

powers.iu <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
powers.or <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
```


####Simulate Results for Politican Response on Support Outcome (Liar's Dividend) -- Study 1 Replication, Full Support Index, Study 2 Estimates
```{r loop of lms to simulate results 7}
#Loop to vary the number of subjects
for (j in 1:length(possible.n)){ N <- possible.n[j] # Pick the sample size

Y0 <- rnorm(n=N, mean=0, sd=1) # control potential outcome 
tau1 <- te_2_iu_study_2 # Hypothesize treatment effect for IU, text
tau2 <- te_2_or_study_2 # Hypothesize treatment effect for OR, text
tau3 <- te_1_iu_video # Hypothesize treatment effect for IU, video
tau4 <- te_1_or_video # Hypothesize treatment effect for OR, video
Y1 <- Y0 + tau1 # treatment potential outcome for IU, text
Y2 <- Y0 + tau2 # treatment potential outcome for OR, text
Y3 <- Y0 + tau3 # treatment potential outcome for IU, video
Y4 <- Y0 + tau4 # treatment potential outcome for OR, video

significant.lm.iu <- rep(NA, sims) # Empty object to count significant experiments 
significant.lm.or <- rep(NA, sims) # Empty object to count significant experiments 

#Inner loop to conduct experiments "sims" times over for each N
for (i in 1:sims){
  
  Z.sim <-  sample(0:4, N, replace = T, prob = c(2/6, 1/6, 1/6, 1/6, 1/6)) # Do a random assignment 
  Y.sim <- case_when(Z.sim == 0 ~ Y0, 
                     Z.sim == 1 ~ Y1,
                     Z.sim == 2 ~ Y2,
                     Z.sim == 3 ~ Y3,
                     Z.sim == 4 ~ Y4) # Reveal outcomes according to assignment 

  fit.lm <- lm_robust(Y.sim ~ factor(Z.sim)) # Do analysis (Simple regression)
  
  p.lm <- summary(fit.lm)$coefficients[2,4] # Extract p-values  
  p.lm2 <- summary(fit.lm)$coefficients[3,4] # Extract p-values 
  significant.lm.iu[i] <- (p.lm <= alpha) # Determine significance according to p <= 0.05
  significant.lm.or[i] <- (p.lm2 <= alpha) # Determine significance according to p <= 0.05
  
}
  powers.iu[j] <- mean(significant.lm.iu) # store average success rate (power) for each N 
  powers.or[j] <- mean(significant.lm.or) # store average success rate (power) for each N 
}

```

####Plot Results -- Study 1 Replication, Full Support Index, Study 2 Estimates
```{r plot results 7}
png("Figures/power_followups_rep_iu_study_2.png")
plot(possible.n, powers.iu, type = "b", main = "Power to Detect Liar's Dividend (Text, IU) by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers.iu, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
#Need over 7000

png("Figures/power_followups_rep_or_study_2.png")
plot(possible.n, powers.or, type = "b", main = "Power to Detect Liar's Dividend (Text, OR) by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers.or, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
#Need 2700
```



#####Set Parameters for Power Analysis
```{r establish parameters for power analysis 8}
alpha <- 0.05 # Standard significance level 
possible.n <- seq(from=2000, to=7000, by=100) # The sample sizes we'll be considering 
sims <- 500 # Number of simulated treatment allocations to conduct for each N 

powers.iu <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
powers.or <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
```


####Simulate Results for Politican Response on Support Outcome (Liar's Dividend) -- Study 1 Replication, Full Support Index, Study 1 + 2 Pooled Estimates
```{r loop of lms to simulate results 8}
#Loop to vary the number of subjects
for (j in 1:length(possible.n)){ N <- possible.n[j] # Pick the sample size

Y0 <- rnorm(n=N, mean=0, sd=1) # control potential outcome 
tau1 <- te_2_iu # Hypothesize treatment effect for IU, text
tau2 <- te_2_or # Hypothesize treatment effect for OR, text
tau3 <- te_1_iu_video # Hypothesize treatment effect for IU, video
tau4 <- te_1_or_video # Hypothesize treatment effect for OR, video
Y1 <- Y0 + tau1 # treatment potential outcome for IU, text
Y2 <- Y0 + tau2 # treatment potential outcome for OR, text
Y3 <- Y0 + tau3 # treatment potential outcome for IU, video
Y4 <- Y0 + tau4 # treatment potential outcome for OR, video

significant.lm.iu <- rep(NA, sims) # Empty object to count significant experiments 
significant.lm.or <- rep(NA, sims) # Empty object to count significant experiments 

#Inner loop to conduct experiments "sims" times over for each N
for (i in 1:sims){
  
  Z.sim <-  sample(0:4, N, replace = T, prob = c(2/6, 1/6, 1/6, 1/6, 1/6)) # Do a random assignment 
  Y.sim <- case_when(Z.sim == 0 ~ Y0, 
                     Z.sim == 1 ~ Y1,
                     Z.sim == 2 ~ Y2,
                     Z.sim == 3 ~ Y3,
                     Z.sim == 4 ~ Y4) # Reveal outcomes according to assignment 

  fit.lm <- lm_robust(Y.sim ~ factor(Z.sim)) # Do analysis (Simple regression)
  
  p.lm <- summary(fit.lm)$coefficients[2,4] # Extract p-values  
  p.lm2 <- summary(fit.lm)$coefficients[3,4] # Extract p-values 
  significant.lm.iu[i] <- (p.lm <= alpha) # Determine significance according to p <= 0.05
  significant.lm.or[i] <- (p.lm2 <= alpha) # Determine significance according to p <= 0.05
  
}
  powers.iu[j] <- mean(significant.lm.iu) # store average success rate (power) for each N 
  powers.or[j] <- mean(significant.lm.or) # store average success rate (power) for each N 
}

```

####Plot Results -- Study 1 Replication, Full Support Index, Study 1 + 2 Pooled Estimates
```{r plot results 8}
png("Figures/power_followups_rep_iu_study_1_2.png")
plot(possible.n, powers.iu, type = "b", main = "Power to Detect Liar's Dividend (Text, IU) by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers.iu, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
#Need 6400

png("Figures/power_followups_rep_or_study_1_2.png")
plot(possible.n, powers.or, type = "b", main = "Power to Detect Liar's Dividend (Text, OR) by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers.or, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
#Need 3100
```



####MDE Calculations for Belief
```{r mde}
#Get SDs of outcome variables using Study 1 data
sds <- ld %>% group_by(alleg_treatment) %>% summarize(belief_sd = sd(belief_1))

n <- seq(from=1000, to=6000, by=250)

belief_MDE <- vector()

for (i in 1:length(n)){
belief_MDE[i] <- 2.8 * sqrt((sds[1,2]^2/(n[i]/3)) + (sds[2,2]^2/(n[i]/3)) + (sds[3,2]^2/(n[i]/3)))
}

#Divide MDEs by control outcome SDs so they are standardized
belief_MDE <- unlist(belief_MDE)/pull(sds[1,2])

png("Figures/mde_followups_belief.png")
mde_data <- as.data.frame(cbind(n, belief_MDE))

mde_plot <- ggplot(mde_data, aes(x=n, y=belief_MDE)) + geom_point() +
  xlab("Sample Size") + ylab("MDE (Standardized)\n") + ggtitle("MDE by Sample Size: Belief 1") +
  # scale_y_continuous(limits = c(.1,.3), breaks = seq(.05, .3, by = .05)) +
  # scale_x_continuous(limits = c(1000,6000), breaks = seq(1000,6000, by = 1000)) +
  theme_bw()
mde_plot

dev.off()
```

